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Abstract 

In this paper, we proposed a new boundary condition scheme for the BGK model 
which has second order accuracy and can be developed to higher accuracy. Numer- 
ical tests show that the numerical solutions of the BGK model applied to the new 
boundary scheme have great agreement with analytical solutions. It is also found 
that the numerical accuracy of the present schemes is much better than that of the 
original extrapolation schemes proposed by Guo et al at small grid. 
Key Words: fluid mechanics; BGK methods; finite difference; boundary conditions; 
numerical method. 
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1 Introduction 

In the past 20 years, the Boltzmann-BGK method in the simulation of fluid showed 
great success. Compared to traditional numerical methods, the Boltzmann-BGK method 
that based on the mesoscopic view got macroscopic variables throughout the moment 
integration of distribution function /. For the mesoscopic property, the Boltzmann-BGK 
can be easily applied to a very large scope, such as multiphase, multicomponent high- 
speed compressible flows and so on[T2l |lH [10]. In simulation, the boundary condition in 
practice is usually given by macroscopic variables, in BGK method we actually use the 
distribution function. Up to now there are no universal method to get the distribution 
function throughout the macroscopic variables, so it's necessary to discuss the distribution 
function of boundary conditions. 

Virtual equilibrium method, interpolation method and non-equilibrium extrapolation 
method are the main boundary schemes used in BGK model nowadays. Most of the 
schemes based on some sort of extrapolation methods to get the distribution functions 
on boundary ^ [T31 El El [2] • In this paper, we got a new boundary scheme based on 
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mathematical analysis. This new boundary scheme has second order accuracy, can be 
easily to develop to higher order accuracy and can be implemented by different viscosity 
between ffow node and the boundary node. In the following sections, we first give implicit- 
explicit scheme of BGK model, and then discuss boundary scheme, in the last, we applied 
this new boundary condition scheme together with the implicit-explicit scheme of BGK 
model to Couette flow and lid driven cavity flow to test accuracy and stability of this 
boundary scheme. 



2 BGK method 

2.1 Basic BGK model 

The standard dynamics theory of mesoscopic model is described by the Boltzmann 
equation as follow 

| + e.V/ = J. (1) 

where / denotes the distribution function, e denotes the molecule velocity vector, and J 
represents the collision term[3]. By the model of BGK, we can easily simplify the collision 
term J to the following formulation 

§|- + e-V/ = ^(/(^'')-/), (2) 

where r denotes the relaxation time and f'^^'^'> denotes the local equilibrium distribution 
function[T]. 

By the Hermite expansion or Taylor expansion, we discrete the infinite velocity space to 
finite space {eo, 62, ■ ■ ■ , e^} to make the equation ([2]) be a set of discrete velocity equations. 

1^ + el ■ V/. = ^(/f^^) - /.), = 0, 1, ■ ■ ■ , N). (3) 

The macroscopic density p and velocity u can be obtained by the moments of the distri- 
bution function 

P = X]/i, pu = ^eifi. (4) 
In 1992, Qian proposed a D2Q9 model|II]. The /^^"'^ of D2Q9 is 



where coi is a weight coefficient, Cg = is the local sound speed. For D2Q9 model, the 

discrete speed is 

/O 1 -1 V2 -V2 -V2 V2\ 

^^ = ^^0 1 -1 72 72 -v/2 -V/2J (^ = 0'l'-'8)' (6) 
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where the c = 1. The sound speed and weight coefficient is 

4/9, |e;|2 = 

=} 



Cs = -^,^^i={ 1/9, |eip = (7) 
1/36, |e;|2 = 2c^ 

Apply Enskog-Chapman expansion to recover the macroscopic incompressible Navier- 
Stokes equations without outer force 

V-M = 0, 

du ^ 1„ (8) 

— + u ■ Vu = — \/p + uAm, 
ot p 

where the /x represents the kinematic viscosity. We need the discrete velocity and equilib- 
rium distribution function satisfying the following equation: 

i i i 

and 

P = rcl (10) 

In this paper, we applied an explicit-implicit scheme to test the validness of the new 
scheme of boundary condition. 

2.2 Time and spatial discrete 

For the calculation of the equations ([s]) of the basic BGK model, we integrate the 
equation ^ in both side to get 

"t+At Of rt+At pt+At 1 

w-l ;(^^"-«- (u) 

/(0,f,e) = /o(x,e). 

where the fo{x,e) represents the initial value. By mean value theorem of integrals, the 
equality ( [ll] ) can be written as 

- /r + el ■ V/.(t^, x)At = ^{ft'\k, x) - Mt^, x)), (12) 

where t < < t+At, the superscript n represents the step number, i.e. /" = f{t, x), f^^^ = 
f{t + At,x). 

For time discrete, we calculates advection term Cj ■ V fi{t^,x)At as an explicit finite- 
difference form and the collision term ^{f'f''\ts_-, x) — fi{t^, x)) as a implicit finite-difference 
form, and we got 
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where the 6 represents the degree of imphcity, in this paper we make it to be 0.5, and 

Z.L Guo and T.S Zhao introduced a new distribution function to remove the imphcity 
of the equation (13), and the new distribution function is 

9^ = f^ + 7r9{f,-ft'^) (14) 



where ir = At/r^j. Applying this new distribution function to equation (13), we can 
delete the imphcity of the collision term and got 



f- 

9-. 
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(eq),n\ 



- -Ate. . V/r + (1 - TT + ne)fj^ + 7r(l - e)ft'^'\ 



(15) 



In the following part, we discuss about the spatial discrete scheme. We apply a mixed- 
difference scheme which combined upwind scheme with central scheme. The central differ- 



ence of IS 

ax a 



The second-order upwind- difference scheme is 



(16) 



1 



2Ax, 



2Ax, 

The mixture form is 



■[^fi{Xa, ■) - 4:fi{Xa - AX„, ■) + fi{Xa - 2AXa, ")] ^/ ^ia ^ 0, 
1 

■[^fi{Xa, ■) - 4/i(x„ + AXa, ■) + fi{Xa + 2AXa, ")] (^ia < 0. 



(17) 



(C 



where ( G [0, 1]. 

By the time discrete and the space discrete, we got the entire iterative form: 



9. 



n+l 



-Ate, ■ V,/r + (1 - TT + vre)/f + 7r(l - 6)^ 
— (^^r + TT^/^^'^)'"), 



{eq),n 



(19) 



where the Vh represents the euqtion (18). By equation (14) and (|4]), we got 



(20) 
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2.3 Scheme of boundary condition 

It's well-known that the scheme of the boundary condition plays a very important role in 
the flow computation. It not only relates to the stability of the computation, but also relates 
to the the accuracy of the computation. In simulation, the boundary condition is usually 
given by macroscopic variables, but in BGK method, we actually use the distribution 
function. For those reasons, we discuss the following boundary conditional scheme for the 
BGK model. Be Inspired by non-equilibrium extrapolation method which was proposed 
by Guo, Zheng and Shi[13j, we designed the new scheme of boundary condition. The basic 
idea of Guo et al is decomposed the distribution function on the boundary into equilibrium 
part and the non-equilibrium part 

Mt,x,) = ft''\t,x,) + ft'\t,x,). (21) 

where the Xp denotes the position of the physic boundary. At present work, we also 
decompose boundary distribution functions into equilibrium part and the non-equilibrium 
part. In BGK model, the equilibrium f-^''\t, Xp) functions can be assumed as the functional 
of u,p,p,T which are related with the time and space [31 For this reason f\^'^\t^Xp) was 
approximated by 

fl''\t,xp) = ft'\t,pit,xf),uit,xp)). (22) 

where the Xf denotes the flow node that next to the Xp i.e xj = Xp + CjAx, and this 
approximation is at least third order accuracy [13]. For the f'l^^'^\t^Xp) part, consider the 
boundary distribution function which satisfies the equation 

^-Mhpil + el . V/.(t, Xp) = ^{ft'\t, Xp) - f,{t, Xp)), (^ = 0, 1, ■ ■ ■ , AT). (23) 

where the r^, is determined by the viscosity between physical boundary and the fluid. If 
we think that the viscosity between physical boundary and the fluid is the same as inner 
flow, we have the r^, = r. In this paper, we set r^, = r. Assume the solution of equation 
(23) can be expanded as 

Mt, Xp) = ft'^it, Xp) + Tj^it, Xp) + r^/f (t, Xp) + ---. (24) 



Substitute (24) into equation (23) and compare the for coefficient r^,, we have 

fl'\t,Xp) = -{^^+e,-V)ft'\t,Xp). (25) 



Based on equation (5| and (22), we can easily got the f'f''\t,Xp) throughout the macro- 
scopic variables. By ff''^\t,Xp), we can got the f'^^\t,Xp) based on equation (25). We 
use T^fi^\t.,Xp) to approximate f^^'^'^\t,Xp) part. By equation (24) , the accuracy of 
f^'^'^\t,Xp) + T^,f-^\t,Xp) to approximate fi{t,Xp) is second order. On the boundary, we 
have distribution 

Mt,xp) = ft'\t,xp) - r,(| + el ■ V)ft'\t,xp). (26) 
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If we want to get the higher accuracy scheme, we can mix Guo's scheme with proposed 
scheme, and we have 



(27) 



+ Mt, Xf) + + el • V)ft'\t, Xf) - ft'\t, Xf), 
To get the value of f'^^^t, Xp), we discrete the time differential operator ^ as 

(^^^^)bt = ^^[ft\t + At, •) - ft'\t, •)], (28) 
and discrete differential operator V at boundary as 

(^^^2^)b5 = ^J^ft'H^a, •) - ^ft'\^a - Ax^, •) + ft'\x^ " 2Ax^, ■)]. (29) 

Summing up the above discussion, we have entire finite difference method for the 
Blotzmann-BGK model. The following part, we will apply above boundary scheme to- 
gether with the implicit-explicit difference scheme to test the validness of the proposed 
boundary scheme. 



3 Numerical Simulations 

In this section, we will apply the new boundary scheme together with the implicit- 
explicit scheme of BGK model to Couette flow and hd driven cavity flow to test accuracy 
and stability of this boundary scheme. 

3.1 Couette flow 

The Couette plate flow is defined in the region 0^x^l,0^?/^l under a periodic 
boundary condition at the entrance and exit. The bottom plate is kept stationary, and 
the top plate moves horizontally with a constant velocity uo- This Couette flow has the 
following analytical solution 

«*(t, y) = (I + 2 J] exp i-i^Xlt) sin (A^y), 0) (30) 

fc=i 

where the = kn/H, k = 1, 2... and H = L = 1.0 . We defined the average error as 



AverageError — 
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where the n denotes the number of the grid. 

We set Re = {Luo)/fi = 10, C = 0.9 and use the grid x A^^ 



5 X 10, X Ny 



40 X 80 for simulation. The proposed scheme 
(|26|) is implemented to top and bottom plates, and the periodic boundary condition is 



10 X 20, X Ny 



20 X 40, N^ X Ny 



implemented to the entrance and exit [T2j . The numerical velocity profiles at time t = 
0.5, t = 5,t = 10,t = 30 together with analytical solutions are plotted in figure [TJ From 
the figure, we can see the numerical solution greatly agreed with the analytical solution. 
The average error of Guo scheme and proposed scheme are plotted in figure [2] Comparing 
Guo's scheme, the proposed scheme at small grid have better accuracy. 




Grid Number, Ny 



Figure 1: Numerical result which applied the Figure 2: Average error of Guo scheme and 

proposed boundary scheme for the couette proposed scheme for couette flow 

flow 



3.2 Two dimensional lid driven square cavity flow 

The lid driven square cavity flow is classical benchmark problem in numerical simulation 
of the fluid. In lid driven square cavity flow, the top lid of the cavity has a constant velocity 
toward right, and the other three boundary hold still. The geometry of the lid driven flow 
is very simply and the phenomena is very complicated. 

In this section we apply proposed boundary scheme for the lid driven cavity flow. The 
width and height is chosen to be L = 1.0. The initial density p = 1.0, top velocity of lid 
is chosen to be m = {ux,Uy) = (0.1,0). The size of the mesh is A^^ x Ny = 128 x 128, set 
At = 0.1 X Ay. The streamline for Re = 400,1000,3000,5000 were plotted in figure [3} 
For Re = 400, we can see three vortices in the figure and the center of the vortex is at 
X = 0.5563, y = 0.6103. For Re = 1000 it also has three vortices and the coordinates of the 
vortex center is x = 0.5334, y = 0.5759. For Re = 3000, we see four vortices in the picture 
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(1)Re =400 



(2) Re = 1000 




(3) Re = 3000 H) = 5000 



Figure 3: Streamline of the lid driven flow at different Reynolds numbers. 
(l)Re=400;(2)Re=1000;(3)Re=3000;(4)Re=5000. 



and the center of the vortex is x = 0.5224, y = 0.5581. When Re = 5000, we can see there 
are four vortices in the flgure and the center of the vortex is x = 0.5210,?/ = 0.5543. All 
four position of the centers of the vortex well matched the paper [121 E]- 

4 Conclusion 

In this paper, we discussed a new boundary scheme together with the implicit-explicit 
scheme for BGK model. The main point of this scheme is to decompose the distribution 
function on the boundary node into equilibrium and non-equilibrium parts. Based on 
the mathematical analysis, we use —T^,ff\t,Xp) to approximate the non-equilibrium. We 
tested this new boundary scheme in Couette flow and lid driven flow. Numerical test 
showed that the numerical solutions of the BGK model applied to proposed boundary 
scheme is very agreement with the former paper [H [l3j and analytical solution and also 
showed second order accuracy. 

The proposed scheme include the different relaxation time that depend on the dif- 
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ferent viscosity between the flow and the boundary. For this reason, we can adjust the 
scheme according to viscosity between boundary and inner flow. By the equation (24), we 
can theoretically develop this scheme to higher order accuracy. The scheme is very simply 
and doesn't need any more node. In one word, comparing with the original extrapolation 
schemes, we have proposed a new method for the implementation of boundary conditions 
for BGK model which are shown to be of second order accuracy, and have better numerical 
stability. 
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